Recoil velocity at 2PN order for spinning black hole binaries 
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We compute the flux of linear momentum carried by gravitational waves emitted from spinning 
binary black holes at 2PN order for generic orbits. In particular we provide explicit expressions of 
three new types of terms, namely next-to-leading order spin-orbit terms at 1.5 PN order, spin-orbit 
tail terms at 2PN order, and spin-spin terms at 2PN order. Restricting ourselves to quasi-circular 
orbits, we integrate the linear-momentum flux over time to obtain the recoil velocity as function 
of orbital frequency. We find that in the so-called superkick configuration the higher-order spin 
corrections can increase the recoil velocity up to a factor ~ 3 with respect to the leading-order PN 
prediction. Whereas the recoil velocity computed in PN theory within the adiabatic approximation 
can accurately describe the early inspiral phase, we find that its fast increase during the late inspiral 
and plunge, and the arbitrariness in determining until when it should be trusted, makes the PN 
predictions for the total recoil not very accurate and robust. Nevertheless, the linear-momentum 
flux at higher PN orders can be employed to build more reliable resummed expressions aimed at 
capturing the non-perturbative effects until merger. Furthermore, we provide expressions valid 
for generic orbits, and accurate at 2PN order, for the energy and angular momentum carried by 
gravitational waves emitted from spinning binary black holes. Specializing to quasi-circular orbits 
we compute the spin-spin terms at 2PN order in the expression for the evolution of the orbital 
frequency and found agreement with Mikoczi, Vasiith and Gergely. We also verified that in the 
limit of extreme mass ratio our expressions for the energy and angular momentum fluxes match 
the ones of Tagoshi, Shibata, Tanaka and Sasaki obtained in the context of black hole perturbation 
theory. 



I. INTRODUCTION 

A. Motivation and summary of results 

In the past few years the study of linear momentum 
carried by gravitational radiation and the subsequent re- 
coil (or kick) velocity it imparts to a binary merger has 
received a lot of attention, as this recoil effect is astro- 
physically very relevant [IJ. There has been a lot of effort 
devoted in quantifying the impact of gravitational recoil 
on stellar mass black hole population and supermassive 
black hole (SMBH) growth scenarios EJ[31ll|5J[i|71IHJ[n], 
as well as on formation of galactic cores [THl E] • In ad- 
dition the observation of a candidate black hole ejected 
from its host galaxy after a merger has recently been re- 
ported [121 . The evidence for a recoiling black hole lies in 
the detection of broad blueshifted emission lines, presum- 
ably from gas carried by the recoiling hole, accompanied 
by a corresponding set of narrow emission lines from the 
gas left behind in the host object. The estimated recoil 
velocity from the line blueshift is ~ 2650 km/s. However 
Refs. [121 E] have proposed alternative scenarios to ex- 
plain the observations of Ref . [12] based on massive black 
hole binary models. 

Recent estimates from binary black hole merger simu- 
lations [TBiisjiniiiHiiiniiinjinjiiiiiiiiiiiisjiisiisT] 

indicate that for some special spin configurations recoil 
velocities of order ~ 4000 km/s could occur in nature. 
Such high kicks are easily strong enough to eject the black 
hole remnant from its host galaxy. Therefore a precise 



understanding of the magnitude of kick velocities and 
their dependence on binary parameters is paramount for 
the development of accurate galactic population synthe- 
sis models and massive black hole formation scenarios. It 
is however currently impossible to simulate the number 
of mergers required to span the expected binary param- 
eter space when spins are included, due to overwhelming 
computational cost. One must therefore develop analyti- 
cal models for the recoil velocity as function of the masses 
and spins of the black holes. 

The first computations of gravitational recoil in binary 
systems were performed by Fitchett ^28j and later by 
Fitchett and Detweiler [33]. These papers relied upon 
earlier work by Peres [30^ and Bekenstein [3l] , who inde- 
pendently computed leading-order expressions for linear 
momentum flux carried by gravitational waves in terms 
of interference between multipole moments of the radia- 
tion field. Fitchett's [28^ result is limited to the regime 
where the binary's dynamics can be accurately described 
by Newtonian physics supplemented by dissipative terms 
due to emission of gravitational waves. It therefore be- 
comes rather inaccurate when the binary is near merger, 
which is where most of the recoil is accumulated. Thus it 
is imperative to include post-Newtonian (PN) corrections 
within this particular framework. To obtain correct recoil 
velocities at higher PN order, one must use the equations 
of motion for the binary at the appropriate PN order, and 
also include additional couplings between multipole mo- 
ments of the radiation field. This extension of the work of 
Fitchett at IPN order for non-spinning binaries has been 
performed by Wiseman [32] . More recently the computa- 
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tion of the recoil for non-spinning binaries at 2PN order 
has been reported by Blanchet, Qusailah and Will f33| 
[henceforth BQW]. 

If the binary contains spinning black holes, then the 
spins of the holes contribute additional terms to the linear 
momentum flux. Kidder |34j has computed the leading- 
order (spin-orbit) contributions from the spins to the re- 
coil. These contributions turn out to be 0.5PN order 
relative to Fitchett's calculation, i.e. to the leading New- 
tonian order, showing that spins, if large, play a cru- 
cial role in determining the recoil, especially since they 
introduce extra asymmetries in the binary. In this pa- 
per we improve the work of Kidder by computing the 
kick velocity including all spin effects up to 2PN order 
beyond Fitchett's leading-order computation. The new 
contributions we compute in this paper are the next- 
to-leading order spin-orbit terms at 1.5PN order, spin- 
orbit tail terms at 2PN order, and the leading spin-spin 
terms at 2PN order. These terms include in particu- 
lar contributions from the quadrupole moment of each 
spinning black hole, which affect both the orbital equa- 
tions of motion at 2PN order, and the time evolution 
of the spins themselves at 1.5PN, through precession in- 
duced by quadrupole-monopole coupling (see for example 
Refs. [35, 36 J. 

The PN computations just described can of course 
claim to provide a reliable estimate of the recoil veloc- 
ity accumulated solely during the early inspiral phase 
preceding the plunge, merger and ringdown. However 
since a significant amount of the total recoil is gener- 
ated during the plunge, merger and ringdown |37l I38j . 
resummation methods and the inclusion of quasi-normal 
modes must be invoked in order to provide a complete 
analytical model of the recoil. Preliminary attempts in 
this direction were pursued in Refs. [371 I3H1 [3S] within 
the effective-one-body approach [401 HH |42] . Our present 
work pushing the calculation of the PN-expanded linear- 
momentum flux at higher orders provides a foundation 
for constructing more accurate resummed versions of the 
linear-momentum flux. 

Our paper is structured as follows. We flrst complete 
our introductory section with a summary of the notation 
and conventions employed throughout. Next in Sec. [IT] we 
provide a somewhat detailed overview of the treatment 
of spins in general relativity and in PN theory. We de- 
fine carefully different spin variables appearing at various 
steps of our computations, e.g. spin variables of the PN 
source multipole moments and spin variables with con- 
stant magnitude. In Sec. |III| we outline the main com- 
putation and give our main results for generic orbits. In 



Sec. IV we specialize our results to quasi-circular orbits 
and integrate the momentum flux to obtain the kick ve- 
locity. We provide numerical estimates of the kick veloc- 
ity accumulated throughout the inspiral for specific con- 
figurations, namely equal mass binaries with spins equal 
in magnitude but opposite in direction. The spins are 
either coUinear with the orbital angular momentum or 
lying in the orbital plane. In the coUinear case we also 



provide an estimate of the kick for equal masses but un- 
equal spins. Finally in Sec.|V] we provide the expressions 
for the fiuxes of energy and angular momentum accurate 
at 2PN order for spinning binary black holes. We provide 
flux expressions for generic orbits and compare with ex- 
isting results in the literature. More specifically we verify 
that in the extreme mass-ratio limit our fluxes match the 
formulas obtained by Tagoshi et al. [43^ in the framework 
of black hole perturbation theory. We also compute the 
spin-spin terms at 2PN order in the expression for the 
evolution of the orbital frequency derived from the usual 
balance argument, and verify that it matches the expres- 
sion obtained by Mikoczi, Vasiith and Gergely 03] when 
one substitutes the proper expression for the quadrupole 
moment of a Kerr black hole, and one neglects contribu- 
tions from magnetic dipoles. 



B. Conventions 

In this paper we consider black holes which can be 
nearly maximally spinning. To reflect this property, our 
PN counting for spin variables is defined as follows. We 
introduce spin variables for each body, say , which are 



related to the true physical spins Sjfue by 



cSi 



(1.1) 



This rescaling stems from the fact that for maximally 
spinning compact objects, the physical spin scales as 
S'truc ~ GM^/c, where M is the body's mass. The con- 
vention for PN order counting thus states that the physi- 
cal spin is of order 0.5PN. On the other hand the rescaled 
spin variables of Eq. (1.1) are of Newtonian order and 



do not contain any hidden power of the speed of light. 
Of course this scaling does not apply to slowly spinning 



objects, for which S^ti 



GM^ 



^spin 



/c . For such bod- 



ies the rescaled spins (1.1) are of 0.5PN order, but such 
slow spins are not targeted by this work. Throughout 
the body of the paper we use geometric units G = c — I. 
However we keep track of the PN order of a given term 
by assigning to it a multiplicative factor which is a power 
of 1/c. This factor should not be thought of as carrying 
dimensions; it appears solely as part of our PN bookkeep- 
ing. We also denote a term of order nPN, i.e. scaling as 
c~^", as being 0{2n). The PN harmonic coordinates are 
denoted as a;^ = {ct, x*), latin indices being spatial. Note 
that the time coordinate a;° carries a PN counting fac- 
tor of c, to track the relative smallness of time variations 
compared to spatial variations in PN theory. In addi- 
tion we denote the antisymmetric permutation symbol 
by {^,zy, A,p}, with {0,1,2,3} = +1. The Lcvi-Civita 
tensor is then 



1 



-.{^i,v,\p}. 



The version with indices down is given by 



(1.2) 



(1.3) 
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which can be derived by simply lowering the indices with 
g^ti/ on the upper-index version. 



II. SPIN VARIABLES IN GENERAL 
RELATIVITY 

In general relativity the covariant treatment of spin is 
somewhat delicate. We propose to start with a discussion 
of a few subtleties that one encounters when dealing with 
spin, in order to hopefully sweep away from the begin- 
ning any potential confusion regarding our analysis, and 
also to provide an intuitive introduction to the topic for 
the reader unfamiliar with these issues. For more infor- 
mation and other recent reviews of treatment of spin in 
general relativity and PN theory, the reader may consult 

Refs. [HlllSlllelllTlllHlligiMllSTllsaisg 



A. Definitions and evolution equations 

Our discussion here relies heavily on the presentation 
by Wald [53], which in turn is based on the works of 
Beiglbock [55|, Madore [56] and Dixon fST]. Consider a 
distribution of matter described by some stress-energy 
tensor T^'^ . At each spacetime event inside the body. 
I.e. any x such that T^"'{x^) ^ 0, one may define a 
spacelike surface 'S{x'^, nP) over the support of T^'' such 
that it is generated by all geodesies orthogonal to a freely 
specifiable timelike unit vector nP defined at x^. [This 
surface is well-defined as long as its generators do not de- 
velop caustics within the matter distribution; we assume 
this is the case for the purpose of our discussion.] We 
next define the total momentum and spin of the matter 
distribution as^ 



that there exists a unique timelike worldline z'^{t) such 



S^^ix^^nP) 



E(a;^,riP) 



(2.1) 

{y-xiPT'^^%y'^)<n:„{y''), 

(2.2) 



where [''■■■'^1 means antisymmetrization with respect to /i 
and V, y'^ are (integration) coordinates on S. It is then 
possible to show that at each event x^ there is a unique 
timelike unit vector qP coUinear with pp, i.e. 



j^Pp''\x^,qP) ^ 0. 



(2.3) 



By identifying nP = qP one selects a preferred spacelike 
surface T,{x^) at each event inside the body. [Henceforth 
we drop the reference to the choice of normal qP in ar- 
guments.] With this result in hand, one can then show 



^ Here we use Riemann normal coordinates at to define the 
integrals. 



that 



p,iz^)SP''{z^) = 0. 



(2.4) 



A set of kinematical constraints on S^'^ like Eq. (2.4) is 
called a spin supplementary condition. The spin sup- 
plementary condition (2.4 1 selects a unique worldline 



adapted to the matter distribution, which is called the 
center-of-mass worldline. Its normalized tangent vec- 
tor is denoted as — dzP/dr. It is the uniqueness 
of the worldline definition (2.4), first suggested by Tul- 
czyjew j58j, that makes it such a natural choice. In the 
literature other choices of spin supplementary conditions 
(and thus different definitions of center-of-mass world- 
lines) are sometimes employed, for example the condition 
u^SP'^ = suggested by Pirani [59] • For a recent review 
of various spin supplementary conditions and their link to 
the selection of a particular center-of-mass worldline^, we 
refer the reader to the paper of Kyrian and Semerak |61j . 
In addition to fixing the worldline, the spin supplemen- 
tary condition (2.4) also reduces the number of indepen- 



dent components of the antisymmetric tensor SP'^ from 
six to three, the correct number of independent spin de- 
grees of freedom expected from Newtonian physics. 

By integrating a Taylor-expanded version of the stress- 
energy conservation equation V^T^"^ = 0, one can 
show [62] that pP and SP" obey 



DpP 

DSP" 
Dt 



= c^[pPu'' ^p'^uP], 



(2.5a) 
(2.5b) 



where Rp^p„ is the Riemann tensor. While the above 
arguments are strictly valid for material bodies with 
non- vanishing stress-energy tensor, it can be shown that 
Eqs. (2.5) can also be applied to black holes (see for ex- 
ample the effective field theory treatment of Porto |46| ) . 
The spin supplementary condition (2.4) also motivates 



the following definition of a spin 1-form 



SP" = 



1 



mc 



(2.6) 



where the (conserved) mass m is defined through pPpf^ — 
—rn^c^. The conservation property of the mass so defined 
is ensured by system (2.5 1. The inverse relation is given 

by 



(2.7) 



^ Intuitively, one can see that a link between the selection of a 
center-of-mass worldline and the definition of the spin of an ob- 
ject must exist, since the split of the total angular momentum 
into "spin" and "orbital" pieces depends explicitly on the choice 
of center-of-mass worldline about which the "spin" piece is de- 
fined, e.g. see Box 5.6 of Ref. I60| . 
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which is obtained by requiring Sf^p^ = 0. Finally we note 
that by taking a r derivative of Eq. (2.4 1, one can derive 



the following relationship between the momentum and 
the 4- velocity 



2c2 



uaXp- 



(2.8) 



By contracting Eq. (2.8 1 with the 4- velocity, one finds 

1/2 



p\u^ = —nic 



Since S^S^"' 



(2.9) 



0, relation (2.8 1 between pf^ and u^^ im- 



plies the exact equivalence between the two conditions 
S^p^ = and S'^w'^ — 0. Thus one can view the re- 
quirement S'^p^ = Sfj_u^ = as stating that in the frame 
instantaneously comoving with the spinning particle, the 
time component of the spin 1-form equals zero. 



B. Post-Newtonian expansion of spin evolution 
equation 

In this paper we specialize to systems where the PN 
expansion is applicable. Our next task is thus the treat- 
ment of system (2.5) in the context of PN theory. For 



our purposes, it is sufficient to analyze the spin evolution 
equation (2.5b I, as the momentum evolution equation 
has already been well studied (see e.g. Ref. |47|) at the 
order we are working at in this paper (2PN). For the 
remainder of this paper we use harmonic coordinates. 

First of all we rewrite Eq. ( 2.5b I into an evolution equa- 
tion for S^. By simply taking a covariant derivative of 
Eq. (2.7 1, it is easy to obtain 



Dt 



27710^""" Dt 
-^,Py,S,^R\p,u-e0''^Pp^S,, (2.10) 



where Eq. (2.5a I has been used. The next step is to ex- 



pand the above evolution equation in the regime where 
PN gravity is valid. Since we are concerned with all 
contributions from spins at 2PN order in the linear mo- 
mentum flux, we need to check if the right-hand side 
of Eq. (2.101, which is quadratic in spin, contributes to 
the precession equations at the order required for our 
computation which is 1.5PN order. To begin with, it is 



clear from Eqs. ( 2.8 1 and ( 2.9 1 that we may replace p^ by 



mcu^ in Eq. (2.101, as the corrections introduced by that 



substitution are well beyond the PN order that interests 
us here. Equation (2.101 then becomes 



Dt 



2^C-M— a/37 

-t- higher PN corrections 



(2.11) 



We next convert proper time derivatives into coordinate 



time derivative to bring Eq. (2.111 closer to the usual 



form of the precession equations. This is accomplished 
using 



Dt 



dT 

c dt 



(2.12) 



where the second line follows from the parametrization 
— w'^(l,i;Yc), u° being determined by normalization 
and = dz^ /dt is the coordinate velocity of the center- 
of-mass worldline 
implies 



Note here that the condition S'^w'^ = 



5o 



-iS*,; 



(2.13) 



which shows that 5o is 0(1). Wc then have 



dSi 
~di 



Si/R Q^pjU 

(2.14) 

As we shall see later in the paper, we require the right- 
hand side of Eq. (2.141 to be accurate to 1.5PN order, or 
0(3). The Ricmann tensor components are at least 0(2) 
(see for example Weinberg [23]), the spatial components 
of the 4-velocity u' are 0(1), u° is 0(0) and Sq is 0(1). 
Therefore the only possibility for the term involving the 
Riemann tensor to contribute in our computation is the 
following: the index v is spatial, the index a is the time 
index, the index A is the time index, and the index p is 
spatial. This implies that the indices /3 and 7 must also 
be spatial due to the antisymmetry of the Levi-Civita 
tensor. Thus the potential contribution comes from the 
components of the Riemann tensor having the structure 



R, 



Ojk- 



However a direct computation shows that these 
components are all 0(3), which, combined with the pres- 
ence of Ui in front, yields a total contribution at 0(4). 
This implies that the leading contributions of the term 
involving the Riemann tensor in the precession equations 
are 0(4) and do not contribute to our computation. 



The remaining term in the right-hand side of Eq. (2.141 



contains the well-known spin-orbit and spin-spin preces 
sion terms. Specializing to a binary system, we obtain 
that the first term on the right-hand side of Eq. (2.141 
gives (for body 1) 



dSi 



(ni2 • ■W12) S'l - 2{vi2 ■ Si)ni2 



TO2 



+{ni2 ■ Si){vi - 2V2) } - 



'3(ni2 • S2)(ni2 X Si) L 



1 



5'2 X Si 



(2.15) 



where V12 — Vi ~ V2 is the relative coordinate velocity, 
ni2 is the unit vector pointing from body 2 to body 1, 
and ri2 is the coordinate orbital separation. The dot 
and cross products are performed with respect to the 
Euclidean spatial metric. The computation presented 
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here assumes the connection coefficients appearing in 
Eq. (2.141 are generated by the other body only. It has 



been shown in Refs. [171 [M] that the divergent self-field 
terms, which arise when one uses delta-function sources 
in the PN field equations, do not contribute to the preces- 
sion equations at the order we are concerned with here. 
This formally justifies our streamlined overview of the 
derivation of the spin evolution equations which simply 
ignores these self-field terms, as one does for example in 
Newtonian physics. Alternatively one could in princi- 
ple avoid using delta-function sources and arrive at the 
same result from a surface integral approach, which re- 
quires knowledge of the vacuum field equations alone (see 
e.g. Ref. [68]). The result (2.151 was originally obtained 
in Ref. [65' as the classical limit of spinning particles in 
quantum field theory. 



The system of equations (2.5) applies to the so-called 
pole-dipole model of an astrophysical object. When 
dealing with systems of Kerr black holes however, sys- 
tem (2.51 is incomplete as one must include, in principle. 



the contributions from all multipole moments of the Kerr 
black holes. However for our computations, we only need 
to include the contribution of the mass quadrupole mo- 
ment to the orbital equations of motion |66l 167] and to 
the precession equations [35j [36]. The precession term 
induced by quadrupole-monopole coupling is |35l 168] 



QM 



3 m2 



(ni2 • 5i)(ni2 X Si). (2.16) 



ties appearing in the precession equations. However, in 
Refs. [Ill [69] the authors computed the "source multi- 
pole moments" in terms of contravariant spin variables. 
For a major part of this work, we stick with the same 
spin variables used in Refs. [371 [HI], and define 



S'^ 1 



2 mB 

c2 ri2 



(2.17) 



Combining Eqs. (2.151, (2.161 and (2.17), we obtain the 



precession equations in the center-of-mass frame in terms 
of the barred spins. They read 



dSi 



7712 



("12 • V12) Si - 2{vi2 ■ Si)ni2 



-{ni2 ■ Si){vi ~ 2V2) 



1 



-3(ni2 • S2){ni2 X Si) 

-3'^{ni2-Si){ni2X Si)l 

mi 



{S2 X Si) 



(2.18) 



In all other sections of the paper, unless otherwise noted, 
the spin variables we use refer to the (contravariant) 
barred spins defined in Eq. (2.17), even though we do 



not carry the bars throughout, for sake of convenience. 
For completeness we also provide the evolution equations 
in the center-of-mass frame, written in terms of the vari- 
ables 



C. Choice of fundamental spin variable 

Based on the discussion of the previous subsection, 
it would seem natural to work with the covariant spin 
variables {SA)i, A = 1,2, as they are the only quanti- 



S1 + S2, 

m ^ 777 - 
J2 Jl- 

7772 '771 



These are 



(2.19a) 
(2.19b) 



— = ^i{vS)-2 — {vA) 

dt c^r^ 777 



3{n- S) + —{n- A) 

771 



„ =, Sm - 
2S+ — A 

771 



3r] 



— f)m — 

4(n-S)+2 — (n- A) 

777 



(n X S) 



2 — (n-5) + (l-47;)(n- A) 
777 



(n X A) 



(2.20a) 



dA 



HTn — — 

2 — (-i;-5) + (-2 + 47/)(-y A) 
777 



-(n-S) + (l-7y)(n- A) 



Sm 



5+(l-27y)A 



A X S' + 3 



2 — {h ■ S) + {1 - Ari){n ■ A) 

777 



(nx S) + 3(1 - 277) 



2{h- S) + —{h- A) 
777 



(n X A)|. 
(2.20b) 



We note that the fundamental spin variables presented literature. Indeed one can check that the spin evolution 
here differ from the ones typically encountered in the 
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equation (2.181 conserves S'^Sf^, but they do not conserve 



as 



the magnitude of a given Kerr hole's spin S'^ S, 
defined in its local asymptotic rest frame. In the litera- 
ture spin variables which conserve the magnitude of the 
Kerr hole's spin are usually preferred, so it is essential to 
relate our spin variables to spins with constant magni- 
tude, which we denote Si 2- In the center-of-mass frame 
the relation between our spin variables and spins with 



constant magnitude was worked out in Refs. [171 [55]. It 
reads 



(_2.21) 

The corresponding transformation rules for S and A are 



A'^ = A 



rjm 



„ ^ Sm - 


rj 




2S+ — A 


V 


m 


2c2 





-S + {l-2ri)A 



1 

2^ 



5m 
1 

m 

dm 



V ■ S + {l-3r])v ■ A 



and the evolution equations for S"^ = + S2 and A'^ = {m/m2)S2 — {m/mi)Si are 



•qm 



dt 2c2r 



2 (n X ti) X 



75^ + 3'^A^ 

m 



377 



A777 

4(n-5^) + 2 — (n- A^ 

m 



(n X S"") 



ntJl 

2 (n-5^) + (l-4?7)(n- A^) 

m 



{h X A'^) ), 



(2.22a) 
(2.22b) 



(2.23a) 



(ri X d) X 



dA'' m 
~dt~ ~ 2c2r2 

+3(1 - 277) 



A777 

3 — S^ + (3-5?7)A^ 
m 



A'^ X 5^= + 3 



2 — (n-5=) + (l-4?7)(n- A^) 

m 



2{n- S') + —{n- A') 

m 



(n X A'^) 



(n X 5"=) 
(2.23b) 



III. LINEAR MOMENTUM FLUX 

The vacuum spacetime surrounding a PN source of 
gravitational radiation can be subdivided into three dis- 
tinct regions [70l [71], each delimited by specific length 
scales. First we have a weak-field near-zone, where 
the PN expansion is valid. The near zone field is 
parametrized by source multipole moments, which encode 
explicit information about the source of the radiation. 
The near-zone extends out to a size < a typical wave- 
length of radiation emitted by the system. At the bound- 
ary of the near-zone begins the local wave zone of the sys- 
tem. The local wave zone is a region of spacetime where 
the effects of background spacetime on wave propagation 
are negligible, i.e. one can describe the gravitational field 
with linearized gravity around a fiat background (asymp- 
totic rest frame of the source) . The gravitational field of 
the local wave zone is parametrized by a set of radiative 
multipole moments, which can be determined in terms 
of the source multipole moments parametrizing the near 
zone. It is during this matching procedure that tail con- 
tributions to the radiative multipole moments of the local 
wave zone can be identified [72 . These tail terms contain 



information on scattering of the waves off the near zone 
curved spacetime. Outside the local wave zone one finds 
the distant wave zone, where one needs to propagate the 
waves on the curved spacetime separating the source and 
the observer. For example when dealing with mergers of 
supermassive black holes at high redshift, cosmological 
effects on wave propagation must be taken into account 
to model the observed waveform. 

In this paper we are concerned with the recoil im- 
parted to the center-of-mass motion of the source due 
to the emission of gravitational waves. Clearly this re- 
coil should be independent of the large scale details of the 
background spacetime into which the source is embedded 
if the size of the source is much smaller than background 
curvature. This is certainly the case for most localized 
astrophysical sources of gravitational waves like binary 
systems, and therefore the physics of the recoil should 
be entirely captured by the interplay between the near 
zone and the local wave zone. The linear momentum 
carried by the waves away from the source is essentially 
due to interference between different radiative multipole 
moments. Thorne |71j gives the complete expression for 
the linear momentum carried by gravitational radiation 
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as an infinite sum of couplings between different radiative 
multipole moments of the local wave zone. 



A. Fundamentals 

At 2PN order the explicit expression of the linear mo- 
mentum flux of Thorne [71 in terms of mass and current 
source multipole moments (/^ = and Jl = Jii...ii 

respectively) is 



dt 



JL7-(4)r(3) J 



16 



r(3) 7(3) 



1 



63 



45 



I - ^ r(4) t(4) j_ jy-i) J 



'kl 



1 



7-(5) j(4) 

1134 



126 



63 



(4) .(3) 
Jjk 



1 r(6) r(5) 



2 r(5) t(5) 

59400 ^^i^^^^i^^"^ 14175 jlmn klmn 



^ r(5) r(4) 



tail terms, 



From this decomposition it should be clear that our com- 
putation requires the orbital equations of motion at 2PN 
order, which include, for spinning binary black holes, 
spin-spin and quadrupole-monopole couplings. Also 
since the third time derivative of Jij is required at 2PN 

1 

accuracy, one needs the time derivative of Jij at 1.5PN 

s 

accuracy, which then implies that the spin precession 
equations are needed at 1.5PN accuracy. To achieve that 
accuracy one needs to include spin-spin and quadrupole- 
monopole couplings in the precession equations in addi- 
tion to the leading-order spin-orbit term [35) 136] . 



The equation of motion in harmonic coordinates for 
the relative orbital separation a; = rn at 2PN order is 
the following [Ml [33] 



(3.1) where 



r^ipN- 



rOso 



ra2PN- 



(3.4) 



where Jr"'' and J)"' denote the vl"^ time derivative of 1l 



and 

and Jl. The tail terms are shown explicitly below in 
Eqs. (3.101 and (3.11). The core of the computation con- 



sists of evaluating the right-hand side of Eq. (3.1) [and 
also the right-hand sides of Eqs. (3.10) and (3.11) be- 



low] for a binary system. More specifically, one needs 
to evaluate the time derivatives of the source multipole 
moments'^, substituting the evolution equations govern- 
ing the binary's dynamics whenever required. The mul- 
tipole moments can be split into orbital contributions 
(non-spinning) and contributions linear in the spins as 
follows** 



CtlPN 



On = — ^w, 



(3.5a) 



(\^■i^)v' ~%^T^^ -2{2 + ^)- 
l r 



(3.5b) 



where 



II 

NS 

II 

s 

Jl 

NS 



Jl 



II = Il + Il, 

NS S 

Jl = Jl + Jl, 

NS S 





II 



1 2 

II 



C NS C NS 



\ Il + 0{5), 

C-^ S 

12 

Jl H — ^ Jl + 

NS C NS 



1 1 
- J 

C s 



1 3 



IL + Oib), 



\ Jl + 0(5), 



Jl + ^ Jl -fO(5). 

s 



(3.2a) 
(3.2b) 

(3.3a) 
(3.3b) 
(3.3c) 
(3.3d) 



-^v X I 7S 
3r 

+ X ( 3S 



Sm 



m 



3 — A 

m 

Sm 



(3.5c) 



a2PN 



m 



?7(3 - 4:r])v^ [v^ --r^ ] - ^r](13 ~ 4:T])—v^ + — r;(l - 3r])r^ - (2 + 25f] + 2T]^)—r^ 
'2/2 r 8 r 



-j(12 + 29^)^ 
4 



n 



77(15 + At])v^ - 3r;(3 + 2T])r^ - (4 + ilf] + 877^) - 



(3.5d) 
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Si ■ S2) - 5{n ■ Si){n ■ S2) 



n+in- S2)Si + ih- Si)S2 }, 



(3.5e) 



2rimr'* 



-Si 



qSl 



-{fl- Sif -bq{h- S2 



-{fl- Si)Si+q{h- S2)S2 



(3.5f) 



Above V = X is the coordinate relative velocity. Equa- 



tion (3.5f ) follows from Ref. [66]. Equations (3.5 1 are only 
valid in the center-of-mass frame of the binary. The ex- 
istence of this center-of-mass frame stems from the fact 
that the general 2PN equations of motion for a binary 
system composed of spinning bodies admit two conserved 
spatial 3-vectors K, and P such that the combination 
Q = K, + Vt can be interpreted as the coordinate loca- 
tion of the center-of-mass. Thus the conserved 3-vector 
T is interpreted as the center-of-mass coordinate velocity. 
It is possible to show [47 that one can perform a (2PN 
accurate) Poincare transformation, under which the 2PN 
equations of motion are invariant, such that in the new 
coordinate systems one has /C'-t-P't = 0. This new coor- 
dinate system is defined as the 2PN center-of-mass frame 
of the binary^ . Note however that if one were to include 
dissipative (radiation-reaction) terms to the equations of 
motion, one would find that "P is not conserved anymore, 
leading to a radiation-reaction induced recoil, which we 
compute here using instead the classic balance argument. 

Given that mi 2 and S\;i are the masses and spins of 
each black hole, the symbols appearing in Eqs. (3.51 are 
defined as 



(3.6a) 
(3.6b) 

(3.6c) 



m = 


■nil + 


bra — 


mi — m2, 




mi 


q = 






m2 ' 




mim2 


77 = 






m^ 


S = 


S1 + S2, 




( S2 


A = 


m 




V™2 



mi 



(3.6d) 
(3.6e) 

(3.6f) 



It is also interesting to note that by introducing the vec- 
tor 



So = 2S 



5m 





h^) Si + 




mi\ 




mi J 







5-2, 
(3.7) 

one may quite neatly combine asiS2 ^^^^ oqm as follows 



ass 



asiS2 ■ 
3 



2mr** 



S^ - 5(n • Sof 



h + 2{n-So)SQ 



(3.8) 



This simple expression can actually be derived easily 
from the spin-spin Hamiltonian for binary black holes 
computed by Damour |36j . which depends solely on the 
spin combination Sq (and orbital elements of course). 
However for objects other than the Kerr black holes of 
general relativity, Eq. (3.8 1 does not hold since the re- 



lationship between their mass quadrupole moment and 
their spin is different from that of a Kerr black hole. 

The tails terms are composed of two main contribu- 
tions, the non-spinning tail terms and the spin-orbit tail 
terms, which we denote as 



dt 



where 



tail 



dt 



NS tail 



dP, 

dt 



, (3.9) 



so tail 



\ dt 



NS tail 



and 



dt 



SO tail 



4m 
63" 
32m 



(5), 



In 



45 



t-T 

2b 
In 



-^^ijk \ (t) I Jki (r) 



In 



11 

+ 12 

t~T 

2b 



t-T 



2b 



dr + lf^{t)i mir) 



In 



dr + 4\t)l lf,\r) 



^31 



t-T 
2b 

In 



In 



97 
60 



2b 



t-T 

2b 



dr 



11 
12 



11 
12 



"1 

(3.10) 



dr 
(3.11) 



9 



In Eq. (3.10), one may use the multipole moments at 



Newtonian order only, i.e. Il ^ II and Jl — > Jl, 

NS NS 

and the time derivatives are evaluated using the New- 
tonian equation of motion. In Eq. (3.111 however, one 



substitutes the 0.5PN expression for the current multi- 

1 

pole moment Jki, i.e. Jki ^ c Jki, and evaluates all 

s 

time derivatives using again the Newtonian equation of 
motion. The spin precession equation is not needed to 
evaluate the spin-orbit tail terms. 



nents along different vectors. For example we write 



dP\ 



dpy 

dt 



dpy 

dt 



V (3.13) 



and provide explicit expressions for each coefficient to 
avoid very lengthy formulas. Again we provide explicit 
results for the instantaneous linear momentum flux only, 
and leave the evaluation of the tail contributions when 
specializing to quasi-circular orbits later on. Our results 
for the instantaneous flux are 



B. Results for generic orbits 

We find the following symbolic structure for the linear 
momentum flux 



dP 

lit 



dp\ 


dt 




1 


' UP 


? 


\~dt 


1 


' UP 


C4 


ydt 



dp 



1 



so 



dP\ 

J IPN 



NLSO 



2PN 



dP 

~dt 



dP\ 

Itt) 



NS tail 



SS 



dP\ 

'^^ J SO tail 



(3.12) 



The new terms that we provide in this paper are the next- 
to-leading spin-orbit terms (NL SO) at 1.5PN, the terms 
quadratic in spins at 2PN order and the spin-orbit tail 
terms at 2PN order. We evaluate the spin-orbit tail terms 
explicitly only when reducing to quasi-circular orbits in 
the next section. 



dP 

~dt 



lOor^ m \ r / 

(3.14a) 



dpy 

N 



I677 m Sm 
dt)^j^ 105r4 ^ 



(25^2 



19f^ 



(3.14b) 

Equations (3.141 match the results of Kidder The 



leading order spin-orbit flux is given by 



dPV 

dt J 



so 



dPy 

dt / 

dp\ 

dt / 



15r5 



IGrj^nv' 2 



ux A 



15r5 

^ — r 

15r5 



[3r{h- A) + 2{v A)], 

(3.15a) 
(3.15b) 

(3.15c) 



The expressions for each flux contribution in Eq. (3.12) Again Eqs. (3.15) match the expression of Kidder [34] . 



can be quite involved, so we split each term into compo- The flrst PN corrections to the flux are 



'dP\ 
dt 



dP 
~dt 



IPN 



IPN 



At] m , Sm 
945r4 m 



6(851 - 77977)w* - 6(2834 - 1877'n)r'^v'^ - 3(4385 - 956r?)— + 6(1843 - 1036?7)f^ 

r 



^(12301 - II6877)— f2 _ g(295 _ 2r;)-^2 



(3.16a) 



Arj m Sm 
945r4 m 



111(25 - 2877)1;'' + 30(392 - 257?7)u^?^-' -I- 9(907 - 162?7)— - 3(2663 - 1394r/)P 

r 



-3(2699 -t- lOr;)— ^ §(^89 + Ut])- 



(3.16b) 



As mentioned in the introduction Wiseman [5T originally cated to compare with our expression, and so we did not 
computed the IPN linear momentum flux, but his results perform that check. We next have the NL SO flux, given 
are presented in a format which makes is quite compli- by 
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'dpy 

. / NLSO 



' 315r5 



226(14 - 53t])v^ - 6(522 - 1973r])r^ + (503 - 1817??)— j [A • (n x v) 



+4^ (l466i;2 - I497f2 + 265^) [S • (n x | , 



(3.17a) 



dPV 
dt 



945r5 



NLSO 



(3968 - 15017?7)t;2 - 3(1274 - 4787r])f'^ + (697 - 2503r?) —J [A • (n x v)] 



+ ^ (7985v^ - 8043f2 + 2176^) [S ■ {h x v)] | 



(3.17b) 



dP 

~dt 
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. 5m 



3f 



4(431 - 2954?7)u2 - 6(348 - 2057r/)r^ + (1508 + 1225r;) — (n • A) 



+12r- 



m 



1628?;^ - 1821r2 _ 212 



TO 



{h-S)+ (-1772 + 48657?)i;2 + 3(1490 - 304l7?)r2 



+ (320- 143r/)-j(t;- A)-2 — 
r J TO 



1531w^ - 2523r^ - 10- 



{v-S)\, 



(3.17c) 



rixA 
'^^ J NL SO 



I 2(661 - 484r?)t;'' - 3(3385 - 2,QMr))r^v^ + 5(233 - 596??) + 9(1143 - 1090r/)f* 



+3(1195+ 150277) I , 



(3.17d) 



4r7 TO (5to / 



1469w* - 7491r^t;^ + 412— i;^ + 6300r* - 1812— - 72- 



,TO 



TO .9 TO 



rff z^NLSO 945r^ TO 



(3.17e) 



/ NL SO 



-^^^f<^ (5677 - 15868r?)w^ - 3(2785 - 6016?7)r^ - 206(29 + 4r/) — 



(3.17f) 



dP_ 
~dt 



vxS 
NL SO 



5431.^- 5709f^ + 1112^ 
945r'' TO r 



(3.17g) 



These next-to-leading order spin-orbit contributions to terms quadratic in spins. We first give the non-spinning 
the Unear momentum flux are new. The contributions at terms 
2PN order contain terms independent of the spins and 
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'dp_ 

dt 



2PN 



10395r4 ~^ 



rl 3(2040 - 18794577 + 149936772)t;6 + 3(42464 + 9003597? - 50304077^)7 



-3(229227 - 45868377 + 17887377^)— _ 24(5363 + 150719?7 - 65604772)7-4^2 

r 

+ (2634273 - 4982252?; + 139140377^)— t-^^^ _,_ (1515394 _ 7543617? + 212216r/ 

r 

+60(84 + 2471377 - 8792772)7-6 _ 3(658810 - 11283917? + 25923677^)— 7-* 



-(1606846 - 56281577 + 866227/2) — /-^ 



2(52781 + 9463877 - 364277^)- 



(3.18a) 



dP 

~dt 



2PN 



277 777 Sm 



18(15482 - 5421577 + 45928772)1;'^ + 18(73439 - 3072407? + 1531807/2)f27;4 



31185r4 m 

+9(141321 - 2148137? + 801737^2)— _ i8(i2l084 - 429935?? + 153642772)7-V 

r 



-9(599979 - 979482?? + 22I85I772)— 7-2^2 _ 2(955835 - 265551?? + 89829?72) — ?;2 

-18(61339 - 177850?? + 48782?;2)?-6 + 3(1448844 - 2083359?? + 336232?/2) —f-'^ 

r 



m 



+6(381131 - 62105?? - 108557?2) —7^2 _ (472594 + 413208?? - 2647277^ 



7?l" 



(3.18b) 



As far as we are aware, expression (3.181 for the 2PN non- 



spinning linear momentum flux for generic orbits has not 
been reported before. When speciahzed to quasi-circular 



orbits, expression (3.181 matches the one of BQW. Fi- 



nally we present the 2PN contributions to the linear mo- 
mentum flux that are quadratic in spins. These terms 
have never been computed before, and are given by 



dPy 

dt 



2^2 



4?? 771 
45^6 



[?-(85 + 987/)(n • A) + (13 - 167?)(t; • A)] 



+2?-(37 + 9877)(n • S) ~ 2(11 + 1677)(i; • S) 
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?-(47 + 19677)(n • A) + 8(7 - 4?7)(^; • A) + 2 f - A9r{n ■ S) 

771 

37;2 - i50f^ + 4—\ [(5 x A) • n] +37-[(S' x A) • i. 



(vS)] 



[{n xv)-S] 
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+2r 



3(463 + 477)7;2 - 3(389 + 887;)?-2 8(5 _ 18??)— {v ■ A) 
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(37- 144??)— (v-S) 

r . 



(3.19b) 



12 



dpy 

dt 



4ri I . Sm 
315r^ 1 m 
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+72f 
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Sm 
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(n-A)^ 
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24(117 - 554r?)7;2 - 3(2471 - 142007?)7-2 - 4(7 + 1447?)- 



(n- A)(n-S') 
(n- A)(t;-S) 



-12r 



Sm 



m 
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m 



(n • S)^ + 12 



Sm 



m 



554^;2 - 1775r2 + 24 



771 



(3.19c) 
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. Sm 



m 
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Sm 
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3(1223 - 44047?)t;2 - 3(2043 - 76167?)f2 + 4(143 - 576??)- 
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(h ■ + 11016f — (n • S)(v ■ S) - 1992 — (t; -Syl, 
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. Sm 



m 



101^2 - 1377^2 
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m 
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Sm 
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- 22f2 
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{vS) 
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'dP\ 
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SS 



45^,6 
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39f^ 



4- [S-{nxv)] + 



9iv^ + 93f^ 



[A • (n X v) 



(3.19f) 



~dt 



tix A 



SS 



15r6 
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38 — A 

m 



(n X t)) . 
(3.19g) 



This concludes the presentation of our results for generic 
orbits. We remind the reader that the spins appearing 
in all the formulas of this section are the (contravariant) 
barred spins of Sec. II C Since the difference between the 
barred spins of Sec. II C| and the more often encountered 
spins with constant magnitude is at IPN order, i.e. an 
0(2) difference, and linear in the spins, only the next- 
lo-leading order spin-orbit linear momentum flux compo- 
nents given by Eqs. (3.171 are affected by this change of 
variables. 



IV. REDUCTION TO APPROXIMATE 
CIRCULAR ORBITS 

In binary systems where spins are dynamically negligi- 
ble, it is well-known that emission of gravitational radi- 
ation pushes the eccentricity of the instantaneous oscu- 
lating orbit toward zero. In PN theory of point-particles 
this osculating orbit is simply found by setting f = and 
solving the resulting equation of motion for the angu- 
lar frequency, which leads to the familiar PN generalized 
Kepler's law. When spins are present however, exact 
circular motion is not a solution to the equations of mo- 
tion generically^. But since the spin-orbit and spin-spin 
accelerations terms responsible for the absence of exact 
circular motion are of 1.5PN and 2PN order respectively, 
it is still expected that the instantaneous osculating or- 
bit of a black hole binary should be nearly circular when 
entering the LIGO band. 

Following Poisson "66], we describe this nearly circu- 
lar motion by treating the spin-dependent acceleration 
terms as a perturbation, and linearize about circular mo- 
tion. This procedure is straightforward and leads to the 
following time-dependent expressions for the orbital sep- 
aration and frequency, which are derived in Appendix [B] 



rit) 



1 



[(A.5S)2-(n.5S)2], (4.1a) 
4^2,2- [(^•^S)'-(--^S)']-(4-lb) 



UJ 



The exception is when spins are coUinear with the orbital angular 
momentum. 



Above r and uj are the orbital averages of r{t) and uj{t), 
and the vector A is given by 



(4.2) 



with Xn — {n X v)/\n x v\, so that h, A and -Ln form 
a right-handed orthonormal basis. In this section we ex- 
press all our results in terms of spin variables with con- 
stant magnitudes |17J IMj- The averages f and w are 
related by a modified version of Kepler's law given by 



(mw)^ 



5— + 3 2 



4to4 



{S^oY - 3(S§ • Xn) 
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Relation (4.3 1 can be inverted to provide the ratio m/r 
as 
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(4.4) 



where x = (mo;)^/'^ following BQW. Lastly the orbital 
velocity is expressed as 



V = f{t)fi + uj{t)r{t) X. 



(4.5) 



At 2PN accuracy we can drop as it is a 4PN quantity, 
i.e. 0(8), and we have 



(4.6) 



Instantaneous linear momentum flux 



To obtain the linear momentum flux in the limit 



of quasi-circular orbits , one substitutes Eqs. (4.11 



Eq. (|4.5[), and then the PN expansion (4.4) into 



Eqs. (pT4p, (3.151, (3.161, (3.171, (3.181 and (3.191. The 
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non-spinning contribution is found to be 



where 
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Pn = 



464772x"/2 



105 



(4.8) 



Equation (4.7 1 matches the instantaneous flux of BQW. 
This provides a good consistency check of our computa- 
tions. The contributions to the hnear momentum flux 
depending on the spins are 
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29 TO 



845 90 , 

1 — "n (A • )(iN • A'^ 

116 29 " ^ 



225 90 (5to /125 45 ^ 

— + (A . A^)(L. ■Sl + -[-^ + ^^,]{X- A=)(Xn • A^) 



„3/2 



(4.9c) 



where we have projected the remaining components along 
S'^, A'^, h X S'^ and h x A'^ on the orthonormal basis 
formed by h, A and L-^. This concludes our discussion 
of the instantaneous linear momentum flux in the limit 
of quasi-circular orbits. 



B. Tail contributions to the linear momentum flux 



The tail contributions to the linear momentum flux 



are formally given by Eqs. (3.101 and (3.111. Since the 



tails contribute at 1.5PN order (non-spinning terms) and 
at 2PN order (spin-orbit terms), the precession dynam- 
ics can be dropped for the purpose of computing these 
terms. One can easily see this as follows. The orbital 
plane precession originates from the spin-orbit acceler- 
ation a so, and therefore introduces 1.5PN relative cor- 
rections to the non-spinning tails. These corrections thus 
contribute at 3PN in the flux. Similarly the spin preces- 
sion dynamics introduce IPN relative corrections to the 
spin-orbit tail, and thus also show up at 3PN in the flux. 
For quasi-circular orbits, the fact that we can ignore pre- 
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cession effects when computing tails implies that we can 
parametrize the unit vectors n and A on a convenient 
time-independent triad as 

n{t) = (cos0(t),sin(/)(t),O), (4.10a) 
\{t) = (-sin(?!)(i),cos(/)(i),0), (4.10b) 

where (t>{t) is the orbital phase as function of time. For 
the purpose of evaluating the tail integrals (3.101 and 



(3.111, it is useful to express the unit vectors fiij) and 



A(r), which depend on the integration variable t, as lin- 
ear combinations of h{t) and \{t). Doing so allows one 
to pull vector quantities outside of the integral over r. 
These linear combinations are 

n(r) = cos[(/)(t) - 0(T)]n(t) - sin[(/)(i) - 0(r)]A(i), 

(4.11a) 

A(t) = sin[0(t) - (/)(r)]n(i) -f'Cos[(/)(i) - ?!)(T)]A(t). 

(4.11b) 

In evaluating the time derivatives of the source multipole 
moments appearing in (3.101 and (3.11 1, it is sufficiently 



accurate to substitute the Newtonian equations of mo- 
tion when necessary. The non-spinning tail contributions 
at 1.5PN order have been reported by BQW, but we re- 
view in detail their computation for sake of completeness, 
and also as a methodology check for our computation of 
the new tail terms involving the spins. We begin with 
the 1.5PN non-spinning tail terms. As a first step, one 
must first compute the index contractions appearing in 
Eq. (3.101. Defining = (p{t) — (/'(r), these are found to 
be 



16 5m f]^ 17/2 

5 TO TO^ 

|202cos(2(/3) A(t) - 203sin(2(y9)n(<)| , 

(4.12a) 

2 Sm jf' 17/0 
—x' X 

5 TO TO 

I [ - cos((/7) -I- 3645 cos(3(/3)] X{t) + 
[sin(95) + 3645sin(3(^)]n(t)| , (4.12b) 

J i 



■' TO TO-^ 

I cos((^) \{t) — sin((/3) n(t)| 



(4.12c) 



I cos(2.^) \{t) + sin(2.^) ri.(i)| . 

(4.12d) 

It is well known [73] that even though the tail integrals 
extend throughout the entire history of the binary, it 



is sufficient here to use the instantaneous Newtonian dy- 
namics of the binary neglecting spin effects and radiation- 
reaction (adiabatic approximation) in order to evaluate 
the tails. Thus we may substitute S(f> = tot — Cut in 



Eqs. (4.121, with the orbital frequency w assumed con- 



stant. To evaluate the t integrals, one only needs the 
following formula 



In f ^) 
\2bJ 



= ^ ( - + i\ln{2Bnuj) + 7e1 V 
nuj L 2 L J J 



(4.13) 

where 7e is the Euler-Mascheroni constant. We provide 
a derivation of this essential expression for the unfamiliar 
reader in Appendix [C] The scale B in the logarithm ker- 
nel for each tail integral appearing in Eq. (3.101 is equal 
to &e~"/i2^ &e-9''/6°, be-'^/^ and be'^^/^^'^pectively. 
Performing the integrations and collecting terms yields 
the result 



'^^ J NS tail 

where 



5m I 3097r . 
TO I ^ 



= Pn—x^''^{ ^^A-f In 



NS 



n >, 



(4.14) 



1 r5921 48, 405, , . x 



Equations (4.14) and (4.151 match the results of BQW. 



We will discuss the terms proportional to the logarithm of 
the orbital frequency in more details shortly. We move on 
to the computation of the tail terms linear in the spins. 
We first provide the index contractions of the relevant 
terms, which are 



,kif{t)j^:hr) 



3^a;^(cos((^) in x A" 

TO* L 



sin(v3) (A X A'^), + [cos((y5) (A • A'^) 
-sin((^)(n. A^)](Ln).}, (4.16a) 

2 

e,jfc4f(t)/jf^(T) = 12^x^1 cos(2^) (n X A^), 

-sin(2(^) (A X A'^), + 
[cos(2(^) (A- A'^) 

+ sin(2(^)(n- A^)](Ln)^}. (4.16b) 

Again using the adiabatic approximation, which here also 
assumes that the spins are kept constant, the errors made 
being of IPN relative order from the spin precession equa- 
tions, the tails integrals can be computed immediately. 
The results are 



In 



37r^(A- A'=)iN + (nx A'^) 
(n- A=)I,N - (A X A'^) 



(4.17) 
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where 



1 (2 

- exp < - ~ 3 In 2 — 7e 



(4.18) 



The leading spin-orbit contribution to the hnear momen- 
tum flux for circular orbits is [see Eqs. (4.9b I, (4.9cl] 



dP\ _p Tx^ 



(A- A^)Ln + («x A^) , (4.19) 



which allows us to rewrite the spin-orbit tail as 



dP\ 

J so tail 



'dP\ 
dt / 



'"SO 



dP\ 

lit J 



so 



(4.20) 



where the derivative c?^ refers to parametrization of the 
vectors h and A in terms of the orbital phase displayed 



in Eqs. (4.101. It becomes clear that the tail terms loga- 
rithmic in frequency can be absorbed in the leading-order 
spin-orbit flux by reparametrizing h and A with a differ- 
ent phase variable ^so deflned as 



SO 



„ mo; , 
2— In 



V^so 



(4.21) 



where we have displayed explicitly the PN scaling of the 
phase modulation. The crucial point is to realize that 
the phase modulation induced by the tail terms is a 4PN 
relative correction to the orbital phase, as one can verify 
by taking a time derivative of Eq. (4.21). Indeed one 
finds 



tpso 



2m 








u ^ 




+ 1 






V^so/ 







(4.22) 



Since u) ~ c~^, the second term above scales as c~®, 
which shows explicitly that it is a 4PN relative correc- 
tion. Since we work only at 2PN order in this paper, 
we ignore this phase modulation henceforth. A similar 
argument is made in BQW regarding the terms that are 
logarithmic in frequency in the non-spinning tails, i.e. 
they can be absorbed into a 4PN phase modulation in 
the leading order non-spinning linear momentum flux. 
This completes our discussion of the tail contributions to 
the linear momentum flux at 2PN order. 



C. Estimate of kick velocity 

We now estimate the kick velocity of spinning black 
hole binaries moving along quasi-circular orbits using the 
2PN linear momentum flux computed in the previous 
section. Since we work at 2PN accuracy, we may con- 
sider the orbital frequency as constant when integrating 
the momentum flux in all terms except the Newtonian 



momentum flux. When integrating the Newtonian mo- 
mentum flux, one must use Eq. (4.1b I for the orbital 
frequency. Next, since the precession equations are of 
the form 0(2) + 0(3), we cannot ignore the time de- 
pendence of the spins in the 0(1) spin-orbit linear mo- 
mentum flux, as that time dependence generates extra 
terms of order 0(3) and 0(4). On the other hand, the 
spins may be considered as constant in the 0(3) and 0(4) 
terms of the linear momentum flux, since the precession 
equations generates extra terms scaling at least as 0(5). 
Similarly we cannot ignore the time dependence of 
due to precession in the 0(1) momentum flux, but it can 
be dropped in the 0(3) and 0(4) pieces. 

To required accuracy, the (indefinite) integrals involv- 
ing fi and A that are needed for computing the kick ve- 
locity are the following 



Xdt 



- 

UJ I 12 



c^2 



(n . SlY - (A . Sg) 



— (n.Ao)(A.io)— ^A 



0(5), 
(4.23a) 



i'^'' dt 



I 



1 

3tD 
2 

3aj 
1 

3a; 



3a; 



Xufe + c'(3), (4.23b) 

X'^h^ + X'^W + X^'^n' 

0(3), (4.23c) 



-n 



ijk 



fi^^X^dt = 



-n 



ijk 



3a; 



1 

3(1' 



(4.23d) 



X'^h'^dt = 



3lu 



1 

3a; 



•X ik \ i I ^ik \ i 

X + n X^ 



3a; 



n*^A'' + 0(3). 



(4.23e) 



Equation (4.23al, required for integrating the Newtonian 
fiux, is obtained by integrating by parts using the exact 
expression A = uj^^{t)fi and using Eq. (4.1bl for uj{t). 



The integrals involving the spin-orbit fiux at 0.5PN 
order are of the form [see Eqs. (4.9b I and (4.9c I ] 



A*A^ L^dt 



X' Al L%dt 



X' Al L'^dt 



0(4). (4.24) 
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One can then substitute the evolution equation for Ln, 

^ . ^1/2 . . 



(4.25) 



and Eqs. (2.23a I and (2.23b I for and A^ respec 



tively, in the integrals (4.241 and perform them explicitly 



by keeping the spins constant and using Eqs. (4.23a I 



(4.23e I. The kick velocity V is obtained from the follow- 



ing expression, 



1 



Pdt. 



(4.26) 



so that V = Vj<iV, we can split the kick velocity into the 
following non-spinning and spin contributions (including 
the tail terms) 



'^s — 



Sm 
m 



1 

71345 



452 



1139 



1] ]x + 



22968 2088 



522 

36761 147101 



309 TT 



.3/2 



68904 



>n, 



(4.28) 



Defining the overall multiplicative factor 

464r/2x4 



105 



where we also included the non-spinning tail term from 
(4.27) Ref. [33], and 
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29 



(4.29b) 
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(4.29c) 



D. Special binary configurations 



We investigate here special mass and spin configura- 
tions for which the recoil velocity has been computed in 
numerical simulations. 



The recoil velocity (4.29) we calculated within the PN 



formalism refers only to that portion of the total recoil 
accumulated during the inspiral phase. As shown in nu- 
merical simulations [TSllIl[I71[Tniiniinil221ll31lll 
l25l [26l [27] , and predicted analytically in Ref. "37] within 
the effective-one-body model [lQl|41], the majority of the 



18 



recoil velocity is produced during the plunge, merger and 
ring-down phases. Quite interestingly, depending on the 
black holes' mass and spin, the integrated recoil velocity 
can reach a peak value (around merger) before decreasing 
to a final, smaller velocity asymptotically. The difference 
between the final kick and the kick at the peak is gen- 
erally denoted as anti-kick. Reference [35] showed that 
the amount of anti-kick depends on the way the different 
modes of the linear momentum flux combine either con- 
structively or destructively during the ring-down phase. 
While Eq. (4.291 only applies to the inspiral portion, if 



anything, pushing Eq. (4.29) until the merger might still 



give a rough estimate of the recoil velocity at the peak, 
which is not necessarily the same as the final, total recoil. 

By contrast, if we are interested in predicting analyti- 
cally and with high accuracy the total recoil velocity we 
cannot rely on the PN-expanded equations (4.291. We 



would need to resum the linear-momentum flux or the 
multipole moments and build non-perturbative expres- 
sions which capture the correct results until merger, and 
augment them by the ringdown phase. This approach is 
followed in the effective-one-body model [371 HOI HI] ■ 



1. Spins colUnear with orbital angular momentum 

If the spins are coUinear with the orbital angular mo- 
mentum, then the projections of the spins along n and 
A vanish, leaving 



V^^ = 0, 
= 0, 



(4.30a) 
(4.30b) 



r.l/2 



29 



[1 + 3^x3/2] (iN • A=) 



470 5m 
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V {Ln- A=) 



10^(Ln • S')' + (11 - 4077) (Xn • S^iL^ • A^) + — f ^ - lOr^ ) {L, 



,3/2 



(4.30c) 



For such configurations the total kick velocity lies entirely 
along n. Let us specialize Eq. (4.301 to an equal mass 
binary Sm — 0, rj = 1/4, for which the individual spins 
are equal in magnitude but opposite in direction, such 



that 5"= 



and lA'^ 



= X = dimensionless spin of 



each individual hole. Since Sm — the non-spinning and 
spin-spin kick contributions vanish, and the total kick 
reduces to 



^ 1 15 35 5 



(4.31) 



where ± denotes whether A*^ is aligned or antialigned 
with Z/N- The recoil velocity for this binary configu- 
ration has been computed in several numerical simula- 
tions [T5J Uni [T71 [23]; since it has neghgible anti-kick 
(e.g., see Fig. 2 in Ref. [IS]), the recoil velocity at the 
peak is close to the final, total recoil. The latter was 
estimated to be ~ ±x450 km/sec. As explained above. 



Eq. (4.31 1 can predict only the recoil velocity around the 
peak. Since the choice of orbital frequency at which to 
compute Eq. (4.311 is rather arbitrary, we choose the 
^ 0.15-0.25, spanning the orbital fre- 



range mcu = x ' 
quencies when in the effective-one-body model the ring- 
down phase is joined to the inspiral phase. The location 
of the latter depends on the black holes' mass and spin. 
(4.31) we obtain |Mcick| 



From Eq 
which brae 



X 114-730 km/sec. 



•cets the numerical result, but has large devi- 
ations from it. 



In Ref. [16, numerical simulations were also carried out 
relaxing the condition S'^ — 0. In this case the magnitude 
for the kick velocity is given by 



IMcickl = l\Xl-X2\x'^'^^Yl 



6 



15 35^ 



29 
420 



(Xi +X2) 



,3/2 



(4.32) 



Notice here that there are non-zero spin-spin contribu- 
tions to the kick. In Ref. [16], it is demonstrated that 
while the measured recoil velocities are not very well ap- 
proximated by the original Kidder kick formula [34' , the 
addition of terms quadratic in spin to the Kidder fitting 
formula provides very good matching to numerical data. 
The fitting function of Ref. [TB] has the following form 



IK 



kick 



\X2 



where 



/(y) = 109.3 - 132.5 y + 23.1 km/s. 



(4.33) 



(4.34) 



Given a fixed X2) the authors of Ref. [T5] show that 
Eqs. (4.33 1-( 4.34| ) help to reduce the fit residuals from 20 
km/s (Kidder formula) to 5 km/s. However it sh ould b e 
pointed out here that the functional form of Eqs. (4.33)- 



(4.34) is not invariant under interchange of particle la- 
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bels''. As was highlighted in Ref. this fundamental 
symmetry provides a guiding principle for building a vi- 
able kick velocity fitting formula over the entire binary 
parameter space. The fit provided in Eq. (4.341 is de- 



rived from a series of simulations where X2 is kept con- 
stant (equal to -0.584) and xi is varied. Thus, there must 
remain some hidden dependence on X2 in the numerical 
coefficients in Eq. (4.341 so that particle-label's symme- 



try is satisfied. Indeed our expression (4.321 suggests the 
following alternative fitting formula 

IVkickl = Ixi - X2|[rfl +C?2(Xl +X2)] 

- \X2\[di{l - y) + X2d2{l - (4.35) 

where y = XilX2 (with \y\ < 1), and where di and d2 
are constants determined from the data of Ref. [16] , the 
results being 



di 
d2 



226.9 km/s, 
67.8 km/s. 



(4.36a) 
(4.36b) 



Note that this modified fit does not change the maximum 
kick value of ~ 450 km/s, obtained when xi = ~X2 — 
±1. Nevertheless Eqs. ( |4.35 )-(4.36 ) may provide better 
results than Eqs. (4.33 1-^^434 1 when both xi and X2 are 
varied. 



2. Spins perpendicular mth orbital angular momentum 

Let us now investigate the so-called superkick config- 
uration studied in several numerical simulations [201 [HI 
m [Ml US (Ml HZ]- We speciahze Eq. ( |4.29| to the case 

0, 7/ = 1/4, for which the 



of an equal mass binary Sm 
individual spins are equal in magnitude but opposite in 
direction, i.e., — Q and \A.'^\/m? — x — dimcnsionlcss 
spin of each individual hole. The spins lie initially on 
the orbital plane. For this particular spin configuration, 
the precession equations (2.23a I and (2.23b I ensure that 



the total spin remains zero and that A*^ remains in the 
orbital plane for all time, precluding precession of the or- 
bital plane. In this case the non-spinning kick vanishes, 
as well as the spin contributions along n and A, leav- 
ing the contribution along L-^ as the lone non-zero term. 
The total kick velocity is thus entirely out of the orbital 
plane and is equal to 



Vkick = -X cos LP x^/'^ { — 
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27r 



.3/2 



2 

15 ■ 120" ■ 5 " J 

(4.37) 

where ip is the angle between A"^ and n. For this bi- 
nary configuration, the anti-kick is absent (e.g., see Fig. 
17 in Ref. ^28^). Thus the recoil velocity at the peak 



^ Note however that Eq. | |4.32| clearly is invariant under a relabel- 
ing of the black holes. 



(around merger) is close to the final, total recoil. As 
done in Sec. jlVD 1 we estimate the recoil velocity at 
the peak from Eq. (4.371 varying the orbital frequency in 
the range mu) ^ 0.15-0.25. Maximizing on ip we obtain 
IVkickl ~ X 357-2300 km/sec. Thus, even the maximum 
value obtained for maximal spins 2300 km/sec is some- 
what below the value 4000 km/s predicted by the numer- 
ical simulation. [Note that due to the fast increase of 
the recoil velocity at high frequency, had we computed 
the recoil at muj ~ .3, w e would have obtained 4527 
km/sec] While Eq. (4.371 for the kick velocity might 



not be quite trustworthy at such high orbital frequen- 
cies, it is worth to note that the higher-order spin terms 
computed in this paper increase the recoil velocity by a 
factor ~ 1.5 — 2.7 with respect to the leading order spin 
term computed by Kidder [31], i.e., with respect to the 
2/15 term in Eq. (4.371. Notice also that the kick pre- 



dicted by PN theory for this superkick configuration is 
linear in the spins, i.e. all spin-spin terms vanish in that 
configuration. 

It is interesting to note that the quasi-circular 
radiation-reaction force in the superkick configuration 
could be deduced from Eqs. (3.151, (3.171 using linear- 



momentum balance arguments. For example at leading 
PN order, Eq. (3.151 says that the radiation-reaction 



force is normal to the orbital plane and changes sign 
as the spins precess on the orbital plane. It reaches its 
maximum value when the spins are coUinear with the in- 
stantaneous orbital velocity vector A, and it is zero when 
the spins are perpendicular to A. This radiation-reaction 
force causes the binary center-of-mass to oscillate with 
increasing amplitude up and down along the direction 
perpendicular to the initial position of the orbital plane. 
The magnitude and direction of the recoil velocity normal 
to the orbital plane is ultimately determined by where in 
the orbit the end of the inspiral (i.e., the merger) occurs. 
This picture was also suggested in Ref. [75] (see Fig. 5 
therein) although there the argument is constructed us- 
ing the conservative dynamics. Because of that reason 
it is not clear to us how the picture of Ref. [75] carries 
over to the radiation-reaction force. However since both 
Ref. [75] and ourselves arrive (qualitatively) the same 
result, we suspect there must exist a relationship be- 
tween the radiation-reaction force driving the recoil and 
the conservative force responsible for frame dragging, al- 
though we are unaware of any explicit formulation of that 
correspondance . 



3. Out-of-plane kick for generic configurations 



Here we rewrite Eq. (4.29c I in terms of individual di 



mensionless spins xa — SA/m\ (omitting the super- 
scripts here for sake clarity of notation below) to provide 
a formula that can be used when comparing to numerical 
simulations, and also to shed some light in the recent con- 
troversy of whether the recoil velocity out-of-plane scales 
like rf or rf 
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If we define the angles G and as follows 

n • A = A-'- cos 6, 
n - S = S-^cos^, 

where 

A^ - |A-(iN-A)LN| 

n2 



(9+1) 



1X2" - Qxil, 



\xi + q\i\, 



(4.38a) 
(4.38b) 



(4.39a) 



(4.39b) 



where Xa ~ Xa — (-^n • Xa)L-m, then the component of 
the kick along the orbital angular momentum axis can be 
rewritten as 



and angular momentum fluxes at 2PN order for spinning 
binary black holes, including all spin contributions. The 
flux formulas, taken from Thorne |71j . are 
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where the tail terms are 



Itt 
dt 



tail 
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2he-V2 
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We note here that throughout this section S and A can 
be freely interchanged with S'^' and A^, as the differences 
generated by this substitution appear only at (relative) 
2.5PN order in the energy and angular momentum fluxes. 



Main results 



The quantities 'Xjy are equal to xa ■ -^n- Thus, the out- 
of-plane recoil velocity has a non trivial dependence on 
T] when higher order PN corrections are included. The 
dominant contribution scales as 77^, but there are addi- 
tional non-negligible contributions scaling as r]^. 



V. ENERGY AND ANGULAR MOMENTUM 
FLUXES 

To provide further checks on our methodology, we pro- 
vide here a (complementary) computation of the energy 



The non-spinning contributions up to 2PN order and 
spin-orbit contributions at 1.5PN order to the energy and 
angular momentum fluxes are all well-known. Our goal 
here is to compute the 2PN terms quadratic in sp ins 
in the energy and angular momentum fluxes, so only a 
handful of terms from Eqs. (5.1) contribute. Since the 



spin-spin acceleration depends solely on the spin combi- 
nation Sq, it turns to be much more natural to write the 
fluxes in terms of S'o and A instead of S and A. ^. The 
results for generic orbits are 
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Test-mass limit 



One important check of our computations is provided 
by the hmiting case of a test mass orbiting a Kerr black 
hole. In Ref. [13], Tagoshi et al. computed a PN ex- 
pansion of the energy flux produced by a test mass in 
a circular equatorial orbit around a Kerr black hole ob- 
tained via the Teukolsky formalism. In this section we 
show that our energy flux matches the expression of [15] 
at 2PN order. Restricting attention to circular orbits in 
the equatorial plane, one can solve for the orbital angu- 
lar frequency uj using Eqs. (3.5) for a non-spinning test 
mass, the result being 
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where x is the dimensionless spin of the Kerr hole, which 
is denoted by q in Ref. [43 . Since it is an observable 
encoded in the gravitational radiation observed at null 
infinity, the orbital frequency is a gauge invariant quan- 
tity. We can therefore use it to relate the harmonic gauge 
radial coordinate r of PN theory to the Boyer-Lindquist 
radial coordinate tq of Ref. j43]. Defining v = {m/roY^^, 
Tagoshi et al. found 
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(5.7) 



The quantity v defined in Ref. |43| is not to bo confused with the 
orbital velocity v of PN theory. 



Equating Eqs. (5.6) and (5.7), one obtains, to 2PN accu 



racy, the following relation 
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Substituting r ~ and v — cor into Eq. (5.4) for the 



energy flux, and supplementing the resulting expression 
with all other contributing terms at 2PN order (see e.g. 
Ref. [48], ignoring however Eq. (F17) for the spin-spin 
energy flux, which is incomplete), one can verify straight- 
forwardly, making use of Eqs. ( |5.7| and (5.8), that the 
resulting energy flux at 2PN order is 



dE 
~dt 



32772 



1 



33 

16' 



1247 . 

V" 

336 

44711 Y 
9072 J 



73 

12-^ 



(5.9) 



which precisely matches Eq. (3.40) of Ref. [43] computed 
from black hole perturbation theory. The 47r term at 
1.5PN order is the contribution from the tail integral 
given by Eq. (5.2). The corresponding computation of 



the angular momentum flux is straightforward, the ad- 
ditional terms contributing at 2PN order being found in 
Refs. [31177]. The result is 
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By using Eq. (5.7) together with Eq. (5.9), one can 



rewrite the angular momentum flux very simply as 
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which verifles Eq. (3.41) of Ref. 
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C. Evolution of the orbital frequency 



where 



We compute here the evolution equation for the orbital 
frequency, derived from the usual energy balance argu- 
ment and specialized to quasi-circular orbits. The bal- 
ance argument says that the (average) orbital frequency 
evolves according to 
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where dE/dt is given by Eq. (S.lal and where £'orb('^) is 



the orbital energy. The orbital energy, which is conserved 
by the 2PN orbital dynamics defined by Eqs. (3.5), is 
given by 
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Reducing Eqs. (5.4| and (5.141 to circular orbits follow- 



ing the same prescription as performed in Sec. IV for 
the linear momentum flux, substituting the results into 
Eq. (5.121 and averaging the result over orbital motion 
yields 
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and where 
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In Eqs. (5.161 and (5.171, we have xa — SA/m^. The 



term involving the sum over A in Eq. (5.17) is generally 



omitted in the literature, but does indeed contribute at 
the same order as the X1X2 piece. It should therefore 
be included in templates for spinning binary black holes. 



Equation (5.17) is equivalent to the sum of Eqs. (9b), (9c) 
and (9d) of Ref. [H], when the quadrupole moment of 
a Kerr black hole is substituted into Eq. (9d). This 
completes our report on energy and angular momentum 
fluxes at 2PN for spinning binaries. 



VI. CONCLUSIONS 

In this paper we computed the linear-momentum flux 
carried by gravitational waves emitted from spinning bi- 
nary black holes at 2PN order for generic orbits, notably 
the next-to- leading order spin-orbit terms at 1.5 PN or- 
der, spin-orbit tail terms at 2PN order, and spin-spin 
terms at 2PN order. In addition, as far as we know, the 
2PN non-spinning terms for generic orbits we provide do 
not seem appear in the literature. We also performed 
the reduction to quasi-circular orbits and integrated the 
simplified flux over time to obtain the kick velocity as 
function of orbital frequency. We specialized our formula 
for the kick velocity of equal mass binary configurations 
where the spins are equal in magnitude, opposite in di- 
rection and are either coUinear with the orbital angular 
momentum or lying in the orbital plane. In particular, 
we found that in the so-called superkick configuration 
the higher-order spin corrections computed in this paper 
can increase the recoil velocity up to a factor ^ 3 with 
respect to the leading-order PN prediction. 

Comparisons between PN and numerical relativity re- 
sults for the gravitational-wave energy fiux have already 
shown that when the latter is computed for quasi-circular 
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orbits in the adiabatic approximation, as done in this 
paper for the linear-momentum flux, it tends to overes- 
timate the numerical energy flux during the last stages 
of inspiral and plunge [78'. Similar conclusions can be 
drawn here for the linear-momentum flux. 

The PN expressions for the linear-momentum flux can 
be used to grasp which asymmetry in the parameter 
space can produce the recoil velocity, or suggest phe- 
nomenological formulas describing numerical-relativity 
results [To t [2l ] l24l [27] . However, not surprisingly, the fast 
increase during the late inspiral and plunge, and the arbi- 
trariness in determining until when those formula should 
be trusted, make the PN predictions not very accurate 
and robust for predicting the recoil velocity at the peak. 
By contrast, the computation of the linear-momentum 
flux at higher PN orders is crucial for building more reli- 
able resummed expressions aimed at capturing the non- 
perturbative effects until merger [37l [39l HQI [41] and pre- 
dict the total recoil velocity. 

We also provided expressions valid for generic orbits, 
and accurate at 2PN order, for the energy and angu- 
lar momentum carried by gravitational waves emitted 
from spinning binary black holes. Specializing to quasi- 
circular orbits we computed the derivative of the orbital 
frequency through 2PN order, and found agreement with 
results of Mikoczi, Vasuth and Gergely [U]. We also 
verified that in the limit of extreme mass ratio our ex- 
pressions for the energy and angular momentum fluxes 
match the results of Tagoshi et al. ^3] obtained in the 
context of black hole perturbation theory. 

It would certainly be quite interesting to extend this 
computation to 3PN order to provide more refined esti- 
mates of the recoil velocity accumulated during the inspi- 



ral and build more accurate resummed expressions. This 
would require the computation of several new source mul- 
tipole moments. In addition the 3PN acceleration and 
3PN spin precession equations currently available in the 
literature [49^ ^ .51] [511 would need to be computed 
in harmonic gauge. 



Acknowledgments 

A.B. and E.R. acknowledge support from NSF grant 
No. PHY-0603762, and A.B. also acknowledges support 
from the Alfred P. Sloan Foundation. L.K. is supported 
by a grant from the Sherman Fairchild Foundation, and 
NSF grants No. PHY-0652952, No. DMS-0553677 and 
No. PHY-0652929. 

The results of this paper were obtained using three 
independent codes based on Mathematica and MathTen- 
sor. 

APPENDIX A: SOURCE MULTIPOLE 
MOMENTS 

Here we list all source multipole moments required for 
our computations. The expressions below are only valid 
in the center-of-mass frame. The mass moments lij and 
lijk, along with the current moment are needed at 
0(4) (2PN) accuracy. The mass moment lijki and the 
current moment Jijk are required at 0(2) (IPN) accu- 
racy, and the moments lijkim and Jijki are needed at 
0(0) (Newtonian) accuracy. We denote with < i. . .j > 
the symmetric trace-free part with respect to the indices 
i and j. The explicit expressions are 
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It is important to note here that in principle there should 
also be a contribution in lij at 2PN order from the indi- 
vidual quadrupole moment of each black hole. Indeed by 
dimensional analysis one finds that the mass quadrupole 
of a Kerr black hole scales as, restoring factors of G and 
c for clarity here, 
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where we used the PN scaling between the true (physical) 
spin and the spin variable with finite limit as c — > cxd dis- 
cussed in the introduction. Since the first time derivative 
of this quadrupole moment comes from the spin preces- 
sion equation only, which is 0(2), the contribution to the 
linear momentum flux from time derivatives of the indi- 
vidual Kerr quadrupoles is pushed to 3PN order, and 
can thus be ignored here. This completes the list of all 
required source multipole moments for the computation 
of the linear momentum flux at 2PN order for spinning 
binaries. 



APPENDIX B: CONSTRUCTION OF 
QUASI-CIRCULAR ORBITS 

When spins are present, exact circular motion in a 
fixed orbital plane is not a solution to equations of mo- 
tion (3.5). The spins induce radial and azimuthal per- 
turbations, as well as precession of the orbital plane. In 



this Appendix we provide, for the benefit of the unfamil- 
iar reader, a review of the well-known fact that despite 
these diSiculties, it is still possible to meaningfully derive 
a modified version of Kepler's law for spinning binaries. 
This modified Kepler's law relates the orbit-averaged or- 
bital frequency and the orbit-averaged orbital separation, 
as given in Eq.(4.3 1. 



Our description of the orbit follows exactly the formal- 
ism of Ref . [79] . The basic picture is the following. One 
describes the orbit using the unit vector h along the line 
joining the two bodies, the unit vector Lj^ normal to the 
instantaneous orbital plane, and the vector A = Z/n x n- 
With respect to this basis, the instantaneous velocity and 
acceleration are shown to be 
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= aN(r, n) + ^aiPN(r', f, n, v) + -ja2p-N{r, r, n, v) 
+ 4aso('', 5', A) + ^ass(?','n., 5'o). (Bib) 

Our quasi-circular orbits are then constructed as follows. 
Since the leading-order spin acceleration is of 1.5PN or- 



der, we assume that the radial perturbations scale simi- 
larly, i.e. r ^ 0(3). Hence at 2PN accuracy we may set 
f = and V = rtoX in the arguments of each accelera- 
tion term in Eq. (Bib I. By projecting the result on the 



(n, A, Xn) triad, we find (setting c = 1) 
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where — r^uP' . Next we decompose r and uj into their 
orbital average piece plus a time-dependent fluctuation, 
i.e. 
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The radial motion and the (time dependent) orbital fre- 
quency are determined from the n and A components of 
the equation of motion (Bib). Let us first look at the 



component of the equation of motion along A. It yields 
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In order to perform the integral, we may keep the spins 
constant, as their time derivatives yield higher order 
terms. We may also use A = ui^^ii in the right-hand 



side of Eq. ( B4 ) , which then gives 
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On the right-hand side we may assume that r and uj 
are constants, as their time derivatives are at least 0(3). 
Hence the 2PN accurate solution to the A component of 
the equations of motion yields 
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where k is an integration constant. Substituting decom- 
position ( B3 1 into Eq. ( B6 1 we find 
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Since, by definition, Sr and Slo have zero orbital average, 
the constant k is determined as 



The 2PN accurate general solution to this difi^erential 
equation is 



K = r uj + 
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where we have used the following orbital average 
(nW) = i(<5^^--i^4) + 0(3). 
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Equation (B7l becomes 
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Let us now look at the n component of the equation of 
motion. It yields 
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where can be taken as constant for the purpose of solv- 
ing Eq. (Bill at 2PN order. Taking the orbital average 
of Eq. (Bill we find, imposing (Sr) = 
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Solving Eq. (B12l for uj at 2PN accuracy, we recover 



Eq. (4.3 1 after performing the replacements r — > r and 



UJ ^ UJ. It remains to solve for the radial perturbation. 
Substituting Eqs. (|B10| and (|B12| into Eq. \Bll\ we 



obtain the following evolution equation for 5r 
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where A and are constants determined by initial condi- 
tions. For definiteness we assume ^ = 0, so that the ho- 
mogeneous solution to the radial perturbation vanishes. 
The last element we need is the angular frequency per- 
turbation, which is given by 
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The perturbations to the orbital frequency and radial 
motion are quadratic in the spins and therefore are 2PN 
corrections. Hence we only need to make the distinction 
between r and f and uj and uj in the 0(0) (Newtonian) 
piece of linear momentum fiux. 



APPENDIX C: COMPUTATION OF TAIL 
INTEGRAL 

We provide here, for completeness and also for the un- 
familiar reader, an explicit derivation of Eq. (4.131, cen- 



tral to the evaluation of tail terms in the limit of quasi- 
circular orbits. 

Consider first the integral 



In X e'*^^ dx 



(CI) 



which can be mapped directly to Eq. (4.131 by a simple 



change of integration variable and a redefinition of fc. To 
evaluate Eq. (CI I, let us examine the following contour 
integral 
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where the contour C is the closed contour obtained by 
the union of the paths Ci, C2, C3 and C4 taken counter- 
clockwise as depicted in Fig. [l] 

On the contour C and inside the region it borders, we 
take In z to be the principal logarithm In z — In r + 
with z = re*^. The principal logarithm is analytic on and 
inside the contour C, and therefore Eq. ( C2 1 vanishes. Let 



us now go through each path which makes up the contour 
of integration. First, the integral from R\ to i?2 on the 
real axis (the path C\ in Fig. [ij matches Eq. (CI) when 
the limits i?i ^ and R2 ^ 00 are taken. [This, of 
course, is the reason why the principal branch of Inz is 
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Integrating by parts the second term in the right hand 
side of Eq. (|C5| yields 
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We next make use of the following integral representation 
for the logarithm of a number approaching zero 



FIG. 1: Contour of integration for evaluating the tail integral. 



chosen to perform the integral (C2).] Next we have the 



contribution from the path C2, i.e., the quarter-circle of 
radius R2 in Fig. [T] This is given by 
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Clearly this integral vanishes in the limit R2 c>o due to 
the presence of the e^^^'^ ^ factor in the integrand. The 
integral of the path C4, i.e., the quarter-circle of radius 
Ri is quite similiar to Eq. ( C3 1 , and is given by 



/ Inze^'^^dz = iRi [ {InR 

JCi J it/ 2 



-kB.i sin 6 ^iO ^0 



(C4) 



Since i?ilni?i ^ as i?i ^ 0, Eq.(|C4| clearly vanishes 
in the limit Ri ^ 0. Lastly the integral over the path C3 
, i.e., the positive imaginary axis, is 



Inze^'^'^dz = 



-Ri 



(lny + j|)e ''y{idy), 



2k 



R2 



-kR2-] 



-I I In ye '^^ dy. 
iRi 



(C5) 



du = In [1 — e "1 ^ In a as a — > 0. 

(C7) 

Thus, taking the limits i?i ^ and i?2 ^ c» we obtain 



fc In y e dy = — \iik - 



00 / ^ — u 



1 -e- 



du. 



(C8) 



The remaining integral can be recognized as the digamma 
function '^{x) = dh\V{x)/dx evaluated at a: = 1. It 
is well known that ^(1) — — 7E) 7e being the Euler- 
Mascheroni constant. Thus, the final result is 



In X e^''^ dx 



r 



In y + i 7 



1 



-ky 



{idy) 



2k 



-Ink + jE 



(C9) 



Setting X — u/2B and k = 2Bnuj in Eq. (C9 1, we recover 
Eq. (|4!l3|. 
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